图1 数据集研发技术路线图
下载摘 要：青藏高原生态脆弱，对气候变化敏感，该区域的土壤侵蚀严重威胁泛第三极生态环境安全。面对目前青藏高原土壤侵蚀（水力侵蚀）评估差异较大的问题，为了反映更加可信的土壤侵蚀空间分布和变化，作者采用融合中国区域地面气象要素驱动数据集、The Climate Hazards group Infrared Precipitation with Stations（CHIRPS）和ERA-Interim的降雨数据、SoilGrids土壤数据、90-m DEM等数据及多种侵蚀因子实现方案（9个R因子×3个K因子×3个LS因子×3个C因子）的集成RUSLE模型对青藏高原1981–2018年土壤侵蚀进行评估，通过"RUSLE-IC-SDR"方法与实测产沙量数据进行比较，形成青藏高原水力侵蚀数据集。该数据集内容包括：（1）集成模型中位数对应的1981–2018年、1981–1998年、1999–2018年多年平均的水力侵蚀数据；（2）1981–2018年、1981–1998年、1999–2018年水力侵蚀变化速率。数据空间分辨率为100 m。数据集储存为.tif格式，由26个数据文件组成，数据量为12.3 GB（压缩为6个文件，5.38 GB）。
关键词：土壤侵蚀；青藏高原；RUSLE集成模型；动态变化；水力侵蚀
DOI: https://doi.org/10.3974/geodp.2025.04.01
CSTR: https://cstr.escience.org.cn/CSTR:20146.14.2025.04.01
数据可用性声明：本文关联实体数据集已在《全球变化数据仓储电子杂志（中英文）》出版，可获取：
https://doi.org/10.3974/geodb.2025.06.05.V1 或 https://cstr.escience.org.cn/CSTR:20146.11.2025.06.05.V1.
青藏高原构造变形活跃，地形及地表物质物理性质复杂，且对气候变化响应敏感[1–3]。这些独特的自然条件和剧烈的气候变化使得该地区的土壤侵蚀正在加剧[3–7]，并形成了大范围、大面积的脆弱地表覆盖层，成为其生态环境建设面临的巨大挑战。剧烈的土壤侵蚀将会对青藏高原的牧业管理、交通运输以及当地居民的生命财产安全构成巨大威胁。基于RUSLE模型的青藏高原土壤侵蚀研究已获得广泛关注。然而，由于所使用数据、侵蚀因子计算方法的不同，不同研究评估的青藏高原水力侵蚀在量级大小、空间分布及动态变化等方面存在显著差异[6,8–11]。这种评估结果的不一致性，不仅阻碍了对区域侵蚀状况形成统一、可靠的认识，也制约了精准水土保持措施的制定和生态环境管理策略的有效实施。因此，亟需对青藏高原水力侵蚀进行精细、稳健、可靠的评估。
RUSLE模型能够充分综合地理、气候等不同的实测资料，而且可以直接利用遥感数据，具有较好的适宜性，是目前使用范围最广的土壤侵蚀模型[12–14]。然而，尽管基于RUSLE模型在青藏高原的土壤侵蚀研究已较为常见，但以往研究多依赖于单一数据源和算法实现各侵蚀因子，其结果存在一定的不确定性和局限性。为克服这一不足，本研究集成多源数据与多算法方案，构建了包含243种土壤侵蚀评估方案（9种R因子×3种K因子×3种LS因子×3种C因子）的集成建模框架。通过融合不同数据源与算法实现的优势，显著提高了土壤侵蚀估算的稳健性与结果可靠性，有效降低因单一数据或单一方法所带来的偏差。这一集成策略增强了对青藏高原水力侵蚀空间分布与动态变化的精细识别能力，生成了1981–2018年土壤侵蚀数据集。该数据结果为区域侵蚀热点识别、政策制定及未来情景预测提供了更为可靠的科学依据[15]。
《集成RUSLE模型评估的青藏高原水力侵蚀数据集（1981–2018）》[16]的名称、作者、地理区域、数据年代、空间分辨率、数据格式、数据量、数据集组成等信息见表1。
| 条目 | 描述 |
|---|---|
| 数据集名称 | 集成RUSLE模型评估的青藏高原水力侵蚀数据集（1981–2018） |
| 数据集短名 | SoilerosionERUSLE_1981-2018 |
| 作者信息 | 黄艳章，中国科学院生态环境研究中心，yzhuang_st@rcees.ac.cn 信忠保，北京林业大学水土保持学院，xinzhongbao@126.com 高光耀，中国科学院生态环境研究中心，gygao@rcees.ac.cn 马英，中国科学院地理科学与资源研究所，maying@igsnrr.ac.cn 杨丽虎，中国科学院地理科学与资源研究所，yanglihu@igsnrr.ac.cn 宋献方，中国科学院地理科学与资源研究所，songxf@igsnrr.ac.cn |
| 地理区域 | 青藏高原 |
| 数据年代 | 1981–2018年 |
| 空间分辨率 | 100 m |
| 数据格式 | .tif |
| 数据量 | 5.38 GB（压缩后） |
| 数据集组成 | 集成模型中位数对应的1981–2018年、1981–1998年、1999–2018年多年平均的水力侵蚀数据；1981–2018年、1981–1998年、1999–2018年水力侵蚀变化速率 |
| 基金项目 | 中华人民共和国科学技术部（2019QZKK0403） |
| 数据计算环境 | Python; ArcGIS |
| 出版与共享服务平台 | 全球变化科学研究数据出版系统 http://www.geodoi.ac.cn |
| 地址 | 北京市朝阳区大屯路甲11号100101，中国科学院地理科学与资源研究所 |
| 数据共享政策 | （1）"数据"以最便利的方式通过互联网系统免费向全社会开放，用户免费浏览、免费下载；（2）最终用户使用"数据"需要按照引用格式在参考文献或适当的位置标注数据来源；（3）增值服务用户或以任何形式散发和传播（包括通过计算机服务器）"数据"的用户需要与《全球变化数据学报（中英文）》编辑部签署书面协议，获得许可；（4）摘取"数据"中的部分记录创作新数据的作者需要遵循10%引用原则，即从本数据集中摘取的数据记录少于新数据集总记录量的10%，同时需要对摘取的数据记录标注数据来源[17] |
| 数据和论文检索系统 | DOI，CSTR，Crossref，DCI，CSCD，CNKI，SciEngine，WDS，GEOSS，PubScholar，CKRSC |
根据RUSLE模型的输入因子，研究需要降雨数据、土壤数据、数字高程（Digital Elevation Model，DEM）数据以及植被数据。研究使用的降雨数据包括3种，分别为中国区域地面气象要素驱动数据集（China Meteorological Forcing Data，CMFD），The Climate Hazards group Infrared Precipitation with Stations数据集（CHIRPS）和ERA-Interim数据[18–20]，数据的时间尺度均为日尺度。土壤数据来源于SoilGrids数据集[21]，本研究所使用的土壤数据包括土壤有机碳、基于美国农业部（United States Department of Agriculture，USAD）的土壤质地、土壤砂粒、黏粒、粉粒、WRB分类等。空间分辨率为90 m的DEM数据来自美国国家航空航天局（National Aeronautics and Space Administration，NASA）航天飞机雷达地形观测任务（SRTM-90 m）[22]。归一化植被指数（Normalized Difference Vegetation Index，NDVI）数据来源于GIMMS-NDVI-3g和MOD13Q1.V6两个数据源[23,24]。本研究为延长植被时间可用性，并保证数据的可比性，采用时空稳定性分析和统计降尺度技术对两个NDVI数据集进行了基于重叠周期（2000年2月–2015年12月）的校正[25]。植被类型数据下载于资源环境科学数据平台的《中国100万植被类型空间分布数据》，来源于《中国植被图集》[26]。研究将所有数据空间分辨率重采样至100 m。
| 编号 | 分类 | 植被类型 | C因子范围 | C因子取值 |
|---|---|---|---|---|
| 1 | 耕地 | 农田 | 0.15–0.2 | 0.175 |
| 2 | 林地 | 针叶林 | 0.0001–0.003 | 0.00155 |
| 3 | 林地 | 混交林 | 0.0001–0.003 | 0.00155 |
| 4 | 林地 | 阔叶林 | 0.0001–0.003 | 0.00155 |
| 5 | 灌木 | 灌木林地 | 0.01–0.15 | 0.08 |
| 6 | 灌木 | 稀疏灌木 | 0.1–0.45 | 0.275 |
| 7 | 草地 | 草原 | 0.05–0.15 | 0.1 |
| 8 | 草地 | 草丛 | 0.05–0.1 | 0.075 |
| 9 | 草地 | 草甸 | 0.01–0.08 | 0.045 |
| 10 | 草地 | 高山植被 | 0.01–0.15 | 0.08 |
| 11 | 其他 | 水体 | 0 | 0 |
| 12 | 其他 | 冰川 | 0 | 0 |
| 13 | 其他 | 裸地 | 0.1–0.55 | 0.325 |
| 14 | 其他 | 湿地 | 0 | 0 |
研究采用在大空间尺度土壤侵蚀评估中应用广泛的RUSLE模型，通过耦合多源数据和多种侵蚀因子计算方法，形成由9种降雨侵蚀力（R）因子、3种土壤可蚀性（K）因子、3种坡长坡度（LS）因子和3种土地覆盖与管理（C）因子相互组合的243种土壤侵蚀评估方案的集成模型。以泥沙连通性（Sediment Connectivity，IC）与泥沙输移比（Sediment Delivery Ratio，SDR）之间的密切关系进一步开展对侵蚀结果的比较，确定基于集成模型的青藏高原土壤侵蚀状况，并以评估结果中值研究青藏高原近40年（1981–2018年）土壤侵蚀动态变化（图1）。
图1 数据集研发技术路线图
下载RUSLE模型是Wischmeier和Smith利用美国30个州长达30年的径流小区资料和人工降雨实验所获得的数据确立的[27]，计算公式如下：
(1)
式中，A是年土壤侵蚀速率（t ha–1 yr–1），R是降雨侵蚀力因子（MJ mm ha–1 h–1 yr–1），K是土壤可蚀性因子（Mg ha h MJ–1 ha–1 mm–1），LS是坡长坡度因子，C是土地覆盖与管理因子，P是土壤保持或预防措施因子。青藏高原由于恶劣的生存环境，人为的土地开发和管理措施相对较少[28]，因此本研究中将P因子忽略（即P设为1）。
R因子的计算基于3个时间分辨率为日的降雨数据集：CMFD、CHIRPS和ERA-interim。同时采用章文波等人（公式2–4）、Xie等人（公式5）和Yin等人（公式6）3种R因子计算方法[29–31]。因此共有9种R因子（即3种降雨数据×3种R因子计算方法）实现方案。
(2)
(3)
(4)
式中，Rj为半月的降雨侵蚀力（MJ mm ha–1 h–1 yr–1）；k为半月时段内的日数（d）；Pd为日雨量（>12 mm）；Pd12为日雨量≥12 mm的日平均雨量（mm）；Py12为日雨量≥12 mm的年平均雨量（mm）。
(5)
式中，Rday为日降雨侵蚀力（MJ mm ha–1 h–1 yr–1）；j为月份（1–12）；Pd为日降雨量（mm）。其中侵蚀和非侵蚀日之间的日降雨量阈值为9.7 mm[31]。
(6)
式中，Rmonth为月降雨侵蚀力（MJ mm ha–1 h–1 yr–1），Pmonth为月降雨量（mm），(Pday)month为该月最大日降雨量（mm）。其中侵蚀和非侵蚀日之间的日降雨量阈值沿用9.7 mm[31]。
K因子的计算目前比较常见的有EPIC（公式7）、Dg（公式8–10）、Nomograph等算法（公式11、12）[27,32,33]，已经得到了很好的应用。根据Tian等人的研究及SoilGrids土壤数据可用性[10]，用于计算土壤可蚀性因子的为土壤表层0–15 cm的数据。
(7)
式中，K是土壤可蚀性因子（Mg ha h MJ–1 ha–1 mm–1）；San指砂粒（0.05–2.0 mm）含量（%）；Sil指粉粒（0.002–0.05 mm）含量（%）；Cla指黏粒（< 0.002 mm）含量（%）；TOC指土壤有机碳含量（%），对于有机质含量在4%以上的土壤样品，采用了4%的上限，以防止对富含有机质土壤可蚀性的低估[27,34]；SN1=1–San/100，0.1317为单位制转换系数。
(8)
(9)
(10)
式中，Dg为几何平均粒径的自然对数；di为第i级最大粒径（mm）；di-1为最小粒径（mm）；fi为相应粒径级的质量分数（%）（如设置黏粒的下限为0.00005 mm，则含有15%黏粒、40%粉粒和45%砂粒土壤的Dg等于–3.569）。
(11)
(12)
式中，N1是粉粒加极细砂含量（0.002–0.1 mm，%），通常取砂粒含量的20%[34]；N2是黏粒含量（<0.002 mm，%）；根据WRB土壤分类，通过查阅表得出土壤结构s；土壤渗透性p根据美国农业部（USDA）的土壤质地分类的土壤渗透性p表查得；M是与土壤质地有关的因素。
本研究LS值的计算采用了自动化地球科学分析系统（System for Automated Geoscientific Analysis，SAGA）软件中的3种算法[15]，分别为Moore、Desmet和Böhner方法[35–37]。首先需要在ArcGIS中对DEM进行填洼处理，在SAGA中使用Flow Accumulation工具计算总流域面积（Total Catchment Area，TCA），使用Flow Width and Specific Catchment Area工具计算特定集水区面积（Specific Catchment Area，SCA），然后使用SAGA的LS factor直接计算LS因子。
C因子是政策和土地使用决策方面最重要的因素，因为它是最容易受人为控制来减少侵蚀的因子[38]。C因子目前的计算方式包括通过卫星遥感数据（如NDVI、植被覆盖率等）计算（公式13）、经验值法以及两种方法的结合（公式14）。
(13)
式中，C为植被因子，NDVI为归一化植被指数，α取值为–2，β取值为1。
某一土地覆盖类型的C因子值是这些土壤流失率的加权平均值，范围在0到1之间。这种经验值最初是通过实地实验得到的，但是通常由于不同的地区、气候等因素往往造成同种植被类型的C因子在某个范围内波动（表2）。本研究第二种C因子计算根据前人的研究和总结，采用经验值的中值来表示某种植被类型的C因子[15]。
第三种C因子计算方法采取Panagos等的方法[38]，结合赋值方法和土地覆盖度进行计算。
(14)
min(Clanduse)是该土地利用类型先前研究的最低值，Range(Clanduse)是指该土地利用类型先前研究所有值的范围，Fcover是植被覆盖度。植被覆盖度越大，C值越小。各植被类型的赋值根据Panagos等和Majhi等的方法[38,39]，并结合青藏高原相关研究综合得到（表2）。
本研究使用R语言中的"cpm"包，应用Kolmogorov-Smirnov检验来识别土壤侵蚀在时间序列中的突变点。突变点累计数量迅速增加的终止年份被用作土壤侵蚀时间序列的突变点，该年份为前期和后期的分界年。与此同时，中国实施的最广泛的植被恢复项目——"退耕还林还草工程"和"退牧还草工程"于1999年启动，这可以作为印证统计方法识别出的突变点的证据[40]。
采用基于像元的线性最小二乘回归计算年土壤侵蚀量的变化趋势[40]。
(15)
式中，Slope是对应变量x时间序列的变化速率；xi是第i年对应的x变量；i是从1到n的整数。
本研究采用"RUSLE-IC-SDR"方法计算了集成模型最小、第一分位数（Q1）、第二分位数（Q2）和第三分位数（Q3）侵蚀结果对应的产沙量。从水利部门（13个）和文献记录（11个）收集的24个水文测量站的观测产沙量数据与集成模型产沙量进行比较，用于最适结果的输出。比较结果（图2）表明，集成模型Q2的评估的泥沙产量与水文站的观测结果更一致（R2=0.81，RMSE=7.04，NSE=0.81）。因此，集成模型的中位数将作为最终输出的数据集。
图2 实测产沙量与集成模型侵蚀结果对应的产沙量的比较图
下载《集成RUSLE模型评估的青藏高原水力侵蚀数据集（1981–2018）》内容包括：（1）集成模型中位数对应的1981–2018年、1981–1998年、1999–2018年多年平均的水力侵蚀数据（单位：t ha–1 yr–1）；（2）1981–2018年、1981–1998年、1999–2018年水力侵蚀变化速率（单位：t ha–1 yr–1 yr–1）。数据空间分辨率为100 m。数据集储存为.tif格式。
青藏高原1981–2018年多年平均土壤侵蚀速率为5.91±2.29 t ha–1 yr–1，即每年有1,526±591 Tg的土壤发生水力侵蚀。在空间上，土壤侵蚀速率较高的区域主要分布在雅鲁藏布江中上游的日喀则、下游的林芝，以及川西、青藏高原东部和帕米尔高原。而柴达木盆地和羌塘高原的土壤侵蚀速率最低（图3）。
图3 1981–2018年青藏高原平均土壤侵蚀速率的空间分布图
（本图基于审图号为GS(2019)1822的标准地图制作，底图无修改）
根据突变点分析和"退耕还林还草"工程实施的时间，作者将1999年作为土壤侵蚀变化的转折点（图4），分为前期（1981–1998年）和后期（1999–2018年）两个时期。青藏高原土壤侵蚀表现为增加趋势（图4a和图5a），后期土壤侵蚀速率显著高于前期（P=9.5e–5），但前期快速增加的趋势有所缓解（前期趋势为0.075 t ha–1 yr–1 yr–1；后期趋势为0.013 t ha–1 yr–1 yr–1）（图4a）。在前期，雅鲁藏布江中游和青藏高原东南部土壤侵蚀趋势增加（图5b）；相反，后期这些地区的土壤侵蚀开始降低（图5c），土壤侵蚀加剧主要集中在青藏高原西部和东部，整个青藏高原土壤侵蚀缓解的面积扩大了65.6%（图5c）。
图4 1981–2018年青藏高原土壤侵蚀速率的时间变化图
下载
图5 青藏高原土壤侵蚀趋势变化图
（本图基于审图号为GS(2019)1822的标准地图制作，底图无修改）
青藏高原是我国重要生态安全屏障，生态环境极其脆弱、成土浅薄，不合理的人类活动和气候变化正在加剧青藏高原水力侵蚀。青藏高原水力侵蚀的研究为环境、生物多样性、自然资源保护提供基础认识。本研究采用集成RUSLE模型，耦合多源数据和多侵蚀因子计算方法，提供了更可靠的土壤侵蚀数据集；识别了土壤侵蚀变化的热点区域，发现土壤侵蚀速率较高的区域主要位于雅鲁藏布江流域和青藏高原东部；并分析了青藏高原1981–2018年土壤侵蚀的变化规律，发现未来青藏高原西部和东部需要重点关注土壤侵蚀的继续加剧。本研究研制的青藏高原水力侵蚀数据可以为针对性的水土保持建设提供支撑和科学指导。
作者分工：宋献方、高光耀、信忠保和黄艳章对数据集的开发做了总体设计；黄艳章和马英设计了模型和算法；宋献方和杨丽虎对论文和数据提出了优化建议；黄艳章撰写了数据论文。
利益冲突声明：本研究不存在研究者以及与公开研究成果有关的利益冲突。